module Tpr

using Dates
using Utils

include("GenQrv.jl")

function load_tpr(tprfilename::AbstractString, settings)
    gmx = settings["environments"]["gmx"]
    wd = dirname(abspath(tprfilename))
    mdp = joinpath(wd, "mdout.mdp")
    err = joinpath(wd, "mdout.err")
    try
        run(pipeline(Cmd(`$gmx dump -quiet -s $tprfilename -nobackup`, dir=wd), stdout=mdp, stderr=err))
    catch
        error_hint("Error: convert tpr failed. Check $err to find out why.")
    end
    return mdp
end

function dimAPBS(file, Imol, minX, maxX, minY, maxY, minZ, maxZ, cfac, fadd, df, gridType, PBEset, PBEset0, PBAset)
    lenX = max(maxX - minX, 0.1)
    cntX = (maxX + minX) / 2
    lenY = max(maxY - minY, 0.1)
    cntY = (maxY + minY)/2
    lenZ = max(maxZ - minZ, 0.1)
    cntZ = (maxZ + minZ)/2
    cX = lenX * cfac
    fX = min(cX, lenX + fadd)
    cY = lenY * cfac
    fY = min(cY, lenY + fadd)
    cZ = lenZ * cfac
    fZ = min(cZ, lenZ + fadd)

    levN=4    # 划分级别
    t=2^(levN+1)
    nX=round(fX/df)-1
    nX=max(t*round(nX/t)+1, 33)
    nY=round(fY/df)-1
    nY=max(t*round(nY/t)+1, 33)
    nZ=round(fZ/df)-1
    nZ=max(t*round(nZ/t)+1, 33)

    if gridType == 0 # GMXPBSA method
        fpre=1
        cfac=1.7
        fX=lenX+2*fadd
        cX=fX*cfac
        nX=t*(Int32(fX/(t*df))+1+fpre)+1
        fY=lenY+2*fadd
        cY=fY*cfac
        nY=t*(Int32(fY/(t*df))+1+fpre)+1
        fZ=lenZ+2*fadd
        cZ=fZ*cfac
        nZ=t*(Int32(fZ/(t*df))+1+fpre)+1
    end

    MGset="mg-auto"
    mem = 200*nX*nY*nZ/1024/1024 # MB

#		npX=nX; npY=nY; npZ=nZ
#		gmem=4000
#		ofrac=0.1
#		if(mem>=gmem) {
#			while(mem>gmem) {
#				maxN=max(npX, max(npY, npZ))
#					 if(maxN==npX) npX = t*((npX-1)/t-1)+1
#				else if(maxN==npY) npY = t*((npY-1)/t-1)+1
#				else if(maxN==npZ) npZ = t*((npZ-1)/t-1)+1
#				mem = 200*npX*npY*npZ/1024./1024
#			}

#			t=nX/npX; if(t>1) npX = int(t*(1+2*ofrac) + 1.0);
#			t=nY/npY; if(t>1) npY = int(t*(1+2*ofrac) + 1.0);
#			t=nZ/npZ; if(t>1) npZ = int(t*(1+2*ofrac) + 1.0);
#			MGset="mg-para\n  ofrac "ofrac"\n  pdime "npX" "npY" "npZ
#		}

    XYZset="  $MGset\n  mol $Imol\n  dime   $nX  $nY  $nZ        # 格点数目, 所需内存: $mem MB" *
        "\n  cglen  $cX  $cY  $cZ        # 粗略格点长度" *
        "\n  fglen  $fX  $fY  $fZ        # 细密格点长度" *
        "\n  fgcent $cntX  $cntY  $cntZ  # 细密格点中心" *
        "\n  cgcent $cntX  $cntY  $cntZ  # 粗略格点中心"

    return "\nELEC name $file\n$XYZset \n$PBEset \nend\n\nELEC name $file~VAC\n$XYZset\n$PBEset0 \nend\n\nAPOLAR name $file~SAS\n  mol $Imol\n$PBAset\nend\n\nprint elecEnergy $file - $file~VAC end\nprint apolEnergy $file~SAS end\n\n"
end

# 根据原子编号获取残基编号
resid(i, res_list, residues) = findfirst(x -> x == res_list[i], residues)

function do_mmpbsa(mdp, trj, ndx, com, pro, lig, ndxfilename, tprfilename, settings)
    # Running settings
    gmx = settings["environments"]["gmx"]
    apbs = settings["environments"]["apbs"]
    MCSH_HOME = devnull

    # mmpbsa parameters
    radType = settings["mmpbsa"]["radType"]
    radLJ0 = settings["mmpbsa"]["radLJ0"]
    meshType = settings["mmpbsa"]["meshType"]
    gridType = settings["mmpbsa"]["gridType"]
    cfac = settings["mmpbsa"]["cfac"]
    fadd = settings["mmpbsa"]["fadd"]
    df = settings["mmpbsa"]["df"]
    dt = settings["mmpbsa"]["dt"]

    # 极性计算设置(Polar)
    PBEset = "
    temp  298.15      # 温度
    pdie  2           # 溶质介电常数
    sdie  78.54       # 溶剂介电常数, 真空1, 水78.54

    npbe              # PB方程求解方法, lpbe(线性), npbe(非线性), smbpe(大小修正)
    bcfl  mdh         # 粗略格点PB方程的边界条件, zero, sdh/mdh(single/multiple Debye-Huckel), focus, map
    srfm  smol        # 构建介质和离子边界的模型, mol(分子表面), smol(平滑分子表面), spl2/4(三次样条/7阶多项式)
    chgm  spl4        # 电荷映射到格点的方法, spl0/2/4, 三线性插值, 立方/四次B样条离散
    swin  0.3         # 立方样条的窗口值, 仅用于 srfm=spl2/4

    srad  1.4         # 溶剂探测半径
    sdens 10          # 表面密度, 每A^2的格点数, (srad=0)或(srfm=spl2/4)时不使用

    ion charge  1 conc 0.15 radius 0.95  # 阳离子的电荷, 浓度, 半径
    ion charge -1 conc 0.15 radius 1.81  # 阴离子

    calcforce  no
    calcenergy comps"

    # 非极性计算设置(Apolar/Non-polar)
    PBAset = "
    temp  298.15 # 温度
    srfm  sacc   # 构建溶剂相关表面或体积的模型
    swin  0.3    # 立方样条窗口(A), 用于定义样条表面

    # SASA
    srad  1.4    # 探测半径(A)
    gamma 1      # 表面张力(kJ/mol-A^2)

    #gamma const 0.0226778 3.84928
    #gamma const 0.027     0
    #gamma const 0.0301248 0         # AMBER-PB4 .0072*cal2J 表面张力, 常数

    press  0     # 压力(kJ/mol-A^3)
    bconc  0     # 溶剂本体密度(A^3)
    sdens 10
    dpos  0.2
    grid  0.1 0.1 0.1

    # SAV
    #srad  1.29      # SAV探测半径(A)
    #press 0.234304  # 压力(kJ/mol-A^3)

    # WCA
    #srad   1.25           # 探测半径(A)
    #sdens  200            # 表面的格点密度(1/A)
    #dpos   0.05           # 表面积导数的计算步长
    #bconc  0.033428       # 溶剂本体密度(A^3)
    #grid   0.45 0.45 0.45 # 算体积分时的格点间距(A)

    calcforce no
    calcenergy total"

    if lig != ""
        withLig=1
    else
        withLig=0
        com = pro
        lig = pro
    end

    pid = "_" * replace(Utils.get_mdparam(mdp, "name")[2:end-1], " " => "_")
    request_hint("Input system name (default: $pid):\n")
    pid = get_input(pid)
    wd = joinpath(dirname(abspath(tprfilename)), "temp")
    request_hint("Input temp dir to place files (default: $wd):\n")
    wd = get_input(wd)
    if !ispath(wd)
        mkdir(wd)
    end
    tprpath = dirname(abspath(tprfilename))
    err = joinpath(wd, "$pid.err")
    qrv = joinpath(wd, "$pid.qrv")
    pdb = joinpath(wd, "$pid.pdb")

    finish_hint("\n0. Finished setting up environmets and parameters.\n\n")

    # 1. 预处理轨迹: 复合物完整化, 团簇化, 居中叠合, 然后生成pdb文件
    trjwho = "$pid~trj_who.xtc"
    trjcnt = "$pid~trj_cnt.xtc"
    trjcls = "$pid~trj_cls.xtc"
    
    # 提取Complex部分
    try
        println("Fixing complex PBC...")
        comcmd = pipeline(`cmd /c echo $com`, 
                          Cmd(`$gmx trjconv -f $trj -s $tprfilename -n $ndxfilename -o $trjwho -pbc whole -dt $dt -nobackup`, dir=tprpath))
        run(pipeline(comcmd, stderr=err))
    catch
        error_hint("Error: gmx trjconv Failed. Check $err to find reasons.\n")
        return
    end

    if withLig == 1
        # useful for single protein and ligand
        println("Centering ligand...")
        ligcmd = pipeline(`cmd /c echo $lig $com`,
                          Cmd(`$gmx trjconv -f $trjwho -s $tprfilename -n $ndxfilename -o $trjcnt -pbc whole -center -dt $dt -nobackup`, dir=tprpath))
        run(pipeline(ligcmd, stderr=err))
        println("Clustering complex trajectory...")
        comcmd = pipeline(`cmd /c echo $com $com`,
                          Cmd(`$gmx trjconv -f $trjcnt -s $tprfilename -n $ndxfilename -o $trjcls -pbc cluster -dt $dt -nobackup`, dir=tprpath))
        run(pipeline(comcmd, stderr=err))
        println("Fitting rotation and translation...")
        ligcmd = pipeline(`cmd /c echo $lig $com`,
                          Cmd(`$gmx trjconv -f $trjcls -s $tprfilename -n $ndxfilename -o $pdb -fit rot+trans -dt $dt -nobackup`, dir=tprpath))
        run(pipeline(ligcmd, stderr=err))
    else
        println("Centering receptor...")
        ligcmd = pipeline(`cmd /c echo $lig $com`,
                          Cmd(`$gmx trjconv -f $trjwho -s $tprfilename -n $ndxfilename -pdb-pbc mol -center -dt $dt -nobackup`, dir=tprpath))
        run(pipeline(ligcmd, stderr=err))
    end
    finish_hint("\n1. Finished pre-processing trajectory.\n\n")

    # 从这里开始, 不再迎合Gromacs的0索引
    pro += 1
    lig += 1
    com += 1

    # 2. 获取每个原子的电荷, 半径, LJ参数, 然后生成qrv文件
    println("Generating qrv file...")
    
    gen_qrv(mdp, ndx, pro, lig; withLig=withLig, qrv=qrv, radType=radType, radLJ0=radLJ0)
    finish_hint("\n2. Finished generating qrv file.\n\n")

    # 3. Mpdb>pqr, 输出apbs, 计算MM, APBS
    println("Running MM-PBSA calculatons...")

    # 获取dt
    # tlist = [parse(Float64, match(r"TITLE.*t=\s*(\d+\.*\d*)\s*step=.*", title).captures[1]) for title in filter(l -> occursin("TITLE", l), readlines(joinpath(pwd(),pdb)))]
    # dt = tlist[2] - tlist[1]

    qrv_file = readlines(qrv)

    # 获取C6和C12参数
    Atyp = parse(Int64, qrv_file[2])
    C6 = zeros(Float64, Atyp, Atyp)
    C12 = zeros(Float64, Atyp, Atyp)
    for i in 1:Atyp
        for j in 1:Atyp
            paralj = split(qrv_file[2 + i])
            C6[i, j] = parse(Float64, paralj[2 * j])
            C12[i, j] = parse(Float64, paralj[2 * j + 1])
        end
    end

    # 获取组信息
    ndxPro = ndx.groups[pro].indexes
    if withLig > 0
        ndxLig = ndx.groups[lig].indexes
    end

    function parseinfo(line)
        line = split(line)
        return [parse(Float64, line[2]) parse(Float64, line[3]) parse(Int64, line[4]) parse(Float64, line[5]) parse(Float64, line[6])]
    end

    # 获取原子类型参数
    mat = [parseinfo(line) for line in qrv_file[Atyp + 3:end]]
    mat = vcat(mat...)
    Qatm = mat[:, 1]
    Ratm = mat[:, 2]
    Catm = mat[:, 3]
    Satm = mat[:, 4]
    Eatm = mat[:, 5]
    resPro = ["P~" * split(line)[end-2] for line in qrv_file[Atyp + 3:end] if split(line)[end] == "Pro"]
    resLig = ["L~" * split(line)[end-2] for line in qrv_file[Atyp + 3:end] if split(line)[end] == "Lig"]

    # 蛋白配体原子数
    Npro = length(resPro)
    Nlig = length(resLig)
    Ncom = Npro + Nlig
    temp = 0
    pdie = 0
    sdie = 0
    Nion = 0
    Qion = zeros(10)
    Cion = zeros(10)

    # 极性参数
    for i in split(PBEset, "\n")
        if i != ""
            paras = split(i)
            if paras[1] == "temp"
                temp = parse(Float64, paras[2])
            end
            if paras[1] == "pdie"
                pdie = parse(Int32, paras[2])
            end
            if paras[1] == "sdie"
                sdie = parse(Float64, paras[2])
            end
            if paras[1] == "ion"
                Nion += 1
                Qion[Nion] = parse(Int32, paras[3])
                Cion[Nion] = parse(Float64, paras[5])
            end
        end
    end

    Iion = sum(Cion .* (Qion .^ 2))
    eps0 = 8.854187812800001e-12
    kb = 1.380649e-23
    Na = 6.02214076e+23
    qe = 1.602176634e-19
    RT2kJ = 8.314462618 * temp / 1E3
    kap = 1E-10 / sqrt(eps0 * kb * temp * sdie / (Iion * qe ^ 2 * Na * 1E3))

    PBEset0 = replace(PBEset, r"sdie.*\d*\.*\d*" => "sdie  1")

    gamma = 0
    _const = 0
    # 非极性参数
    for i in split(PBAset, "\n")
        if i != ""
            paras = split(i)
            if endswith(paras[1], "gamma") && paras[2] == "const"
                global gamma = parse(Float64, paras[3])
                global _const = parse(Float64, paras[4])
            end
        end
    end

    # 处理TITLE行
    Nfrm = 0
    curpdblines = readlines(pdb)
    titleindex = findall(l -> startswith(l, "TITLE"), curpdblines)
    timestamp = [parse(Float64, match(r".*t=\s*(\d+\.*\d*)\s*.*", curpdblines[t]).captures[1]) / 1000 for t in titleindex]
    Fname = ["$pid" * @sprintf("~%.0fns", t) for t in timestamp]

    minXpro = zeros(length(titleindex))
    minXlig = zeros(length(titleindex))
    minYpro = zeros(length(titleindex))
    minYlig = zeros(length(titleindex))
    minZpro = zeros(length(titleindex))
    minZlig = zeros(length(titleindex))

    maxXpro = zeros(length(titleindex))
    maxXlig = zeros(length(titleindex))
    maxYpro = zeros(length(titleindex))
    maxYlig = zeros(length(titleindex))
    maxZpro = zeros(length(titleindex))
    maxZlig = zeros(length(titleindex))

    minXcom = zeros(length(titleindex))
    minYcom = zeros(length(titleindex))
    minZcom = zeros(length(titleindex))
    maxXcom = zeros(length(titleindex))
    maxYcom = zeros(length(titleindex))
    maxZcom = zeros(length(titleindex))

    # 处理ATOM行

    function parse_atom(line)
        return [strip(line[1:6]), parse(Int64, strip(line[7:11])), line[13:16], strip(line[18:20]), strip(line[22:22]) != "" ? line[22:22] : "X", parse(Int64, line[23:26]), parse(Float64, line[31:38]), parse(Float64, line[39:46]), parse(Float64, line[47:54])]
    end

    Nfrm = length(titleindex)

    minX = 0
    maxX = 0
    minY = 0
    maxY = 0
    minZ = 0
    maxZ = 0

    for i in 1:Nfrm
        pqr_com = open(joinpath(wd, Fname[i] * "_com.pqr"), "w")
        pqr_pro = open(joinpath(wd, Fname[i] * "_pro.pqr"), "w")
        pqr_lig = open(joinpath(wd, Fname[i] * "_lig.pqr"), "w")

        atomlines = filter(l -> startswith(l, "ATOM"), curpdblines[titleindex[i]:(i + 1 > length(titleindex) ? end : titleindex[i + 1])])
        atoms = [parse_atom(atom) for atom in atomlines]

        minXpro[i] = min([atom[7] for atom in atoms if atom[2] in ndxPro]...)
        maxXpro[i] = max([atom[7] for atom in atoms if atom[2] in ndxPro]...)
        minYpro[i] = min([atom[8] for atom in atoms if atom[2] in ndxPro]...)
        maxYpro[i] = max([atom[8] for atom in atoms if atom[2] in ndxPro]...)
        minZpro[i] = min([atom[9] for atom in atoms if atom[2] in ndxPro]...)
        maxZpro[i] = max([atom[9] for atom in atoms if atom[2] in ndxPro]...)
        if withLig == 1
            minXlig[i] = min([atom[7] for atom in atoms if atom[2] in ndxLig]...)
            maxXlig[i] = max([atom[7] for atom in atoms if atom[2] in ndxLig]...)
            minYlig[i] = min([atom[8] for atom in atoms if atom[2] in ndxLig]...)
            maxYlig[i] = max([atom[8] for atom in atoms if atom[2] in ndxLig]...)
            minZlig[i] = min([atom[9] for atom in atoms if atom[2] in ndxLig]...)
            maxZlig[i] = max([atom[9] for atom in atoms if atom[2] in ndxLig]...)
        end
        
        for atom in atoms
            ATOM, INDX, NAME, RES, CHN, NUM, X, Y, Z = atom
            r = Ratm[INDX]
            atomline_pqr = @sprintf("%-6s%5d %-4s %3s %s%4d    %8.3f %8.3f %8.3f %12.6f %12.6f\n", ATOM, INDX, NAME, RES, CHN, NUM, X, Y, Z, Qatm[INDX], r)
            if INDX in ndxPro
                write(pqr_pro, atomline_pqr)
                minXpro[i] = min(minXpro[i], X - r)
                maxXpro[i] = max(maxXpro[i], X + r)
                minYpro[i] = min(minYpro[i], Y - r)
                maxYpro[i] = max(maxYpro[i], Y + r)
                minZpro[i] = min(minZpro[i], Z - r)
                maxZpro[i] = max(maxZpro[i], Z + r)
            end

            if withLig == 1
                write(pqr_com, atomline_pqr)
                if INDX in ndxLig
                    write(pqr_lig, atomline_pqr)
                    minXlig[i] = min(minXlig[i], X - r)
                    maxXlig[i] = max(maxXlig[i], X + r)
                    minYlig[i] = min(minYlig[i], Y - r)
                    maxYlig[i] = max(maxYlig[i], Y + r)
                    minZlig[i] = min(minZlig[i], Z - r)
                    maxZlig[i] = max(maxZlig[i], Z + r)
                end
            end
            
            minXcom[i] = min(minXpro[i], minXlig[i])
            maxXcom[i] = max(maxXpro[i], maxXlig[i])
            minYcom[i] = min(minYpro[i], minYlig[i])
            maxYcom[i] = max(maxYpro[i], maxYlig[i])
            minZcom[i] = min(minZpro[i], minZlig[i])
            maxZcom[i] = max(maxZpro[i], maxZlig[i])
        
        end

        minX = min(minXcom...)
        maxX = max(maxXcom...)
        minY = min(minYcom...)
        maxY = max(maxYcom...)
        minZ = min(minZcom...)
        maxZ = max(maxZcom...)

        close(pqr_com)
        close(pqr_pro)
        if withLig == 1
            close(pqr_lig)
        end
    end

    kJcou = 1389.35457520287
    Rcut = 1E10              # large enough

    residues = [sort([Set(resPro)...])..., sort([Set(resLig)...])...]
    Nres = length(residues)
    
    resPBSA = open(joinpath(tprpath, "$pid~resPBSA.dat"), "w")
    resPBSA_PB = open(joinpath(tprpath, "$pid~resPBSA_PB.dat"), "w")
    resPBSA_SA = open(joinpath(tprpath, "$pid~resPBSA_SA.dat"), "w")
    resMM_VDW = open(joinpath(tprpath, "$pid~resMM_VDW.dat"), "w")
    resMM = open(joinpath(tprpath, "$pid~resMM.dat"), "w")
    resMM_COU = open(joinpath(tprpath, "$pid~resMM_COU.dat"), "w")
    res_MMPBSA = open(joinpath(tprpath, "$pid~res_MMPBSA.dat"), "w")
    MMPBSA = open(joinpath(tprpath, "$pid~MMPBSA.dat"), "w")

    head = "#Frame   " * join([@sprintf("%10s", replace(r, r"~0+" => "~")) for r in residues]) * "\n"
    headDH = "#Frame   " * join([@sprintf("%21s", replace(r, r"~0+" => "~") * "(with DH)") for r in residues]) * "\n"
    write(resPBSA, head)
    write(resPBSA_PB, head)
    write(resPBSA_SA, head)
    if withLig == 1
        write(resMM_VDW, head)
        write(resMM, headDH)
        write(resMM_COU, headDH)
        write(res_MMPBSA, headDH)
    end

    write(MMPBSA, "   #Frame      Binding( with DH ) |    MM    ( with DH )    PB        SA     |   COU    ( with DH )     VDW   |       PBcom        PBpro        PBlig  |    SAcom     SApro     SAlig\n")

    maxstr = max([length(fname) for fname in Fname]...)

    dE = zeros(Float64, Nres)
    dGres = zeros(Float64, Nres, 2)
    dHres = zeros(Float64, Nres, 2)
    MMres = zeros(Float64, Nres, 2)
    COUres = zeros(Float64, Nres, 2)
    VDWres = zeros(Float64, Nres)
    dPBres = zeros(Float64, Nres)
    dSAres = zeros(Float64, Nres)

    vdw = zeros(Nfrm)
    pb = zeros(Nfrm)
    sa = zeros(Nfrm)
    cou = zeros(Nfrm, 2)
    mm = zeros(Nfrm, 2)
    dh = zeros(Nfrm, 2)

    PBres = zeros(Nres)
    SAres = zeros(Nres)

    preK = -1
    if withLig == 1
        preK = 1
    end
    Ipro = ndx.groups[pro].indexes[1]
    Ilig = ndx.groups[lig].indexes[1]

    type = []
    res = []
    x = []
    y = []
    z = []
    resID = []

    for fr in 1:Nfrm
        t_start = now()
        flush(stdout)

        Fout = Fname[fr]
        println(@sprintf("Running frame %3d/%d: %-10s", fr, Nfrm, Fout))

        inputpqr = joinpath(wd, Fout * "_pro.pqr")
        if withLig == 1
            inputpqr = joinpath(wd, Fout * "_com.pqr")
        end
        type = [split(line)[3] for line in readlines(inputpqr)]
        res = [split(line)[4] for line in readlines(inputpqr)]
        x = [parse(Float64, split(line)[end - 4]) for line in readlines(inputpqr)]
        y = [parse(Float64, split(line)[end - 3]) for line in readlines(inputpqr)]
        z = [parse(Float64, split(line)[end - 2]) for line in readlines(inputpqr)]
        resID = [parse(Int32, replace(string(split(line)[end - 5]), r"[A-Z]+" => "")) for line in readlines(inputpqr)]

        # MM
        dEcou = zeros(Float64, Nres)
        dEcouDH = zeros(Float64, Nres)
        dEvdw = zeros(Float64, Nres)
        if withLig == 1
            for i in 1:Npro
                ii = i + Ipro - 1
                qi = Qatm[ii]
                ci = Catm[ii]
                si = Satm[ii]
                ei = Eatm[ii]
                xi = x[ii]
                yi = y[ii]
                zi = z[ii]
                for j in 1:Nlig
                    jj = j + Ilig - 1
                    cj = Catm[jj]
                    r = sqrt((xi - x[jj]) ^ 2 + (yi - y[jj]) ^ 2 + (zi - z[jj]) ^ 2)
                    if r < Rcut
                        t = 1 / (0.1 * r) ^ 6
                        Ecou = qi * Qatm[jj] / r
                        EcouDH = Ecou * exp(-kap * r)
                        Evdw = (C12[Int64(ci) + 1, Int64(cj) + 1] * t - C6[Int64(ci) + 1, Int64(cj) + 1]) * t
                        dEcou[findfirst(x -> x == resPro[i], residues)] += Ecou
                        dEcou[findfirst(x -> x == resLig[j], residues)] += Ecou
                        dEcouDH[findfirst(x -> x == resPro[i], residues)] += EcouDH
                        dEcouDH[findfirst(x -> x == resLig[j], residues)] += EcouDH
                        dEvdw[findfirst(x -> x == resPro[i], residues)] += Evdw
                        dEvdw[findfirst(x -> x == resLig[j], residues)] += Evdw
                    end
                end
            end

            for i in eachindex(dEcou)
                dEcou[i] *= kJcou / (2 * pdie)
                dEcouDH[i] *= kJcou / (2 * pdie)
                dEvdw[i] /= 2
            end
        end
        Evdw = sum(dEvdw)
        Ecou = sum(dEcou)
        EcouDH = sum(dEcouDH)

        # PBSA
        f = open(joinpath(wd, "$Fout.apbs"), "w")
        if withLig == 1
            write(f, "read\n")
            write(f, "  mol pqr $(Fout)_com.pqr\n")
            write(f, "  mol pqr $(Fout)_pro.pqr\n")
            write(f, "  mol pqr $(Fout)_lig.pqr\n")
            write(f, "end\n\n")
        else
            write(f, "read\n")
            write(f, "  mol pqr $(Fout)_pro.pqr\n")
            write(f, "end\n\n")
        end

        if meshType == 0 # GMXPBSA
            if withLig == 1
                write(f, dimAPBS("$(Fout)_com", 1, minX, maxX, minY, maxY, minZ, maxZ, cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
                write(f, dimAPBS("$(Fout)_pro", 2, minX, maxX, minY, maxY, minZ, maxZ, cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
                write(f, dimAPBS("$(Fout)_lig", 3, minX, maxX, minY, maxY, minZ, maxZ, cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
            else
                write(f, dimAPBS("$(Fout)_pro", 1, minX, maxX, minY, maxY, minZ, maxZ, cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
            end
        elseif meshType == 1 # g_mmpbsa
            if withLig == 1
                write(f, dimAPBS("$(Fout)_com", 1, minXcom[fr], maxXcom[fr], minYcom[fr], maxYcom[fr], minZcom[fr], maxZcom[fr], cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
                write(f, dimAPBS("$(Fout)_pro", 2, minXpro[fr], maxXpro[fr], minYpro[fr], maxYpro[fr], minZpro[fr], maxZpro[fr], cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
                write(f, dimAPBS("$(Fout)_lig", 3, minXlig[fr], maxXlig[fr], minYlig[fr], maxYlig[fr], minZlig[fr], maxZlig[fr], cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
            else
                write(f, dimAPBS("$(Fout)_pro", 1, minXpro[fr], maxXpro[fr], minYpro[fr], maxYpro[fr], minZpro[fr], maxZpro[fr], cfac, fadd, df, gridType, PBEset, PBEset0, PBAset))
            end
        end
        close(f)

        try
            run(pipeline(Cmd(`$apbs $Fout.apbs`, dir=wd), stdout=joinpath(wd, "$Fout.out"), stderr=joinpath(wd, "$Fout.out")))
        catch
            error_hint("Error running apbs. Check $(joinpath(wd, "$Fout.out")) for details.\n")
        end

        apbs_out = readlines(joinpath(wd, "$Fout.out"))
        Esol = zeros(3, Ncom)
        Evac = zeros(3, Ncom)
        Esas = zeros(3, Ncom)
        ln = 1
        while ln <= length(apbs_out)
            line = apbs_out[ln]
            if contains(line, "CALCULATION #")
                if contains(line, "($(Fout)_com")
                    t = 1
                    n = Ncom
                end
                if contains(line, "($(Fout)_pro")
                    t = 2
                    n = Npro
                end
                if contains(line, "($(Fout)_lig")
                    t = 3
                    n = Nlig
                end
                if contains(line, "~VAC)")
                    t += 10
                end
                if contains(line, "~SAS)")
                    t += 20
                end
                while true
                    ln += 1
                    line = apbs_out[ln]
                    if (t < 20 && contains(line, "Per-atom energies:")) || (t > 20 && contains(line, "Solvent Accessible Surface Area"))
                        break
                    end
                end

                for j in range(1, n)
                    ln += 1
                    line = apbs_out[ln]
                    paras = split(line)
                    if t < 20
                        r = parse(Float64, paras[3])
                    else
                        r = parse(Float64, paras[end])
                    end
                    if t < 10
                        Esol[t % 10, j] = r
                    elseif t < 20
                        Evac[t % 10, j] = r
                    elseif t < 30
                        Esas[t % 10, j] = gamma * r + _const / n
                    end
                end
            end
            ln += 1
        end
        
        PBcom = 0
        SAcom = 0
        PBpro = 0
        SApro = 0
        PBlig = 0
        SAlig = 0
        for i = 1:Ncom
            Esol[1, i] -= Evac[1, i]
            PBcom += Esol[1, i]
            SAcom += Esas[1, i]
        end
        for i = 1:Npro
            Esol[2, i] -= Evac[2, i]
            PBpro += Esol[2, i]
            SApro += Esas[2, i]
        end
        for i = 1:Nlig
            Esol[3, i] -= Evac[3, i]
            PBlig += Esol[3, i]
            SAlig += Esas[3, i]
        end

        for i in range(1, Npro)
            dPBres[resid(i, resPro, residues)] += Esol[1, Ipro + i - 1] - Esol[2, i]
            dSAres[resid(i, resPro, residues)] += Esas[1, Ipro + i - 1] - Esas[2, i]
        end
        for i in range(1, Nlig)
            dPBres[resid(i, resLig, residues)] += Esol[1, Ilig + i - 1] - Esol[3, i]
            dSAres[resid(i, resLig, residues)] += Esas[1, Ilig + i - 1] - Esas[3, i]
        end

        vdw[fr] = Evdw
        pb[fr] = preK * (PBcom - PBpro - PBlig)
        sa[fr] = preK * (SAcom - SApro - SAlig)
        cou[fr, 1] = Ecou
        cou[fr, 2] = EcouDH
        mm[fr, 1] = cou[fr, 1] + vdw[fr]
        mm[fr, 2] = cou[fr, 2] + vdw[fr]
        dh[fr, 1] = preK * mm[fr, 1] + pb[fr] + sa[fr]
        dh[fr, 2] = preK * mm[fr, 2] + pb[fr] + sa[fr]

        write(MMPBSA, @sprintf("%-12s %9.3f(%9.3f) | %9.3f(%9.3f) %9.3f %9.3f | %9.3f(%9.3f) %9.3f | %12.3f %12.3f %12.3f | %9.3f %9.3f %9.3f\n",
        Fout, dh[fr, 1], dh[fr, 2], mm[fr, 1], mm[fr, 2], pb[fr], sa[fr], 
        Ecou, EcouDH, Evdw, PBcom, PBpro, PBlig, SAcom, SApro, SAlig))

        for i = 1:Nres
            ii = ""
            if i == 1
                ii = @sprintf("%-9s", Fout)
            end
            endres = ""
            if i == Nres
                endres = "\n"
            end
            if withLig == 1
                write(resMM_VDW, @sprintf("%s%9.3f %s", ii, dEvdw[i], endres))
                write(resMM_COU, @sprintf("%s%9.3f(%9.3f) %s", ii, dEcou[i], dEcouDH[i], endres))
                write(resMM, @sprintf("%s%9.3f(%9.3f) %s", ii, dEcou[i] + dEvdw[i], dEcouDH[i] + dEvdw[i], endres))
                write(res_MMPBSA, @sprintf("%s%9.3f(%9.3f) %s", ii, dEcou[i] + dEvdw[i] + dPBres[i] + dSAres[i], dEcouDH[i] + dEvdw[i] + dPBres[i] + dSAres[i], endres))

                COUres[i, 1] += dEcou[i]
                COUres[i, 2] += dEcouDH[i]
                VDWres[i] += dEvdw[i]
            end
            PBres[i] += dPBres[i]
            SAres[i] += dSAres[i]
            write(resPBSA, @sprintf("%s%9.3f %s", ii, preK * (dPBres[i] + dSAres[i]), endres))
            write(resPBSA_PB, @sprintf("%s%9.3f %s", ii, preK * (dPBres[i]), endres))
            write(resPBSA_SA, @sprintf("%s%9.3f %s", ii, preK * (dSAres[i]), endres))
        end

        resPBSA_pdb = open(joinpath(wd, "$pid~resPBSA.pdb"), "w")
        resPBSA_PB_pdb = open(joinpath(wd, "$pid~resPBSA_PB.pdb"), "w")
        resPBSA_SA_pdb = open(joinpath(wd, "$pid~resPBSA_SA.pdb"), "w")
        resMM_VDW_pdb = open(joinpath(wd, "$pid~resMM_VDW.pdb"), "w")
        resMM_pdb = open(joinpath(wd, "$pid~resMM.pdb"), "w")
        resMM_COU_pdb = open(joinpath(wd, "$pid~resMM_COU.pdb"), "w")
        res_MMPBSA_pdb = open(joinpath(wd, "$pid~res_MMPBSA.pdb"), "w")

        for i = 1:Npro
            ii = Ipro + i - 1
            atomline = @sprintf("ATOM  %5d %-4s %3s A%4d    %8.3f%8.3f%8.3f", ii, type[ii], res[ii], resID[ii], x[ii], y[ii], z[ii])

            write(resPBSA_PB_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, preK * dPBres[resid(i, resPro, residues)], preK * dPBres[resid(i, resPro, residues)]))
            write(resPBSA_SA_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, preK * dSAres[resid(i, resPro, residues)], preK * dSAres[resid(i, resPro, residues)]))
            if withLig == 1
                write(resMM_COU_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEcou[resid(i, resPro, residues)], dEcouDH[resid(i, resPro, residues)]))
                write(resMM_VDW_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEvdw[resid(i, resPro, residues)], dEvdw[resid(i, resPro, residues)]))
                write(resMM_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEcou[resid(i, resPro, residues)] + dEvdw[resid(i, resPro, residues)], dEcouDH[resid(i, resPro, residues)] + dEvdw[resid(i, resPro, residues)]))
                write(resPBSA_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dPBres[resid(i, resPro, residues)] + dSAres[resid(i, resPro, residues)], dPBres[resid(i, resPro, residues)] + dSAres[resid(i, resPro, residues)]))
                write(res_MMPBSA_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEcou[resid(i, resPro, residues)] + dEvdw[resid(i, resPro, residues)] + preK * (dPBres[resid(i, resPro, residues)] + dSAres[resid(i, resPro, residues)]), 
                dEcouDH[resid(i, resPro, residues)] + dEvdw[resid(i, resPro, residues)] + preK * (dPBres[resid(i, resPro, residues)] + dSAres[resid(i, resPro, residues)])))
            else
                write(resPBSA_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, preK * (dPBres[resid(i, resPro, residues)] + dSAres[resid(i, resPro, residues)]), preK * (dPBres[resid(i, resPro, residues)] + dSAres[resid(i, resPro, residues)])))
            end
        end
        for i = 1:Nlig
            ii = Ilig + i - 1
            atomline = @sprintf("ATOM  %5d %-4s %3s A%4d    %8.3f%8.3f%8.3f", 
                    ii, type[ii], res[ii], resID[ii], x[ii], y[ii], z[ii])
            write(res_MMPBSA_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEcou[resid(i, resLig, residues)] + dEvdw[resid(i, resLig, residues)] + preK * (dPBres[resid(i, resLig, residues)] + dSAres[resid(i, resLig, residues)]), 
                dEcouDH[resid(i, resLig, residues)] + dEvdw[resid(i, resLig, residues)] + preK * (dPBres[resid(i, resLig, residues)] + dSAres[resid(i, resLig, residues)])))
            write(resPBSA_PB_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, preK * dPBres[resid(i, resLig, residues)], preK * dPBres[resid(i, resLig, residues)]))
            write(resPBSA_SA_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, preK * dSAres[resid(i, resLig, residues)], preK * dSAres[resid(i, resLig, residues)]))
            write(resMM_COU_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEcou[resid(i, resLig, residues)], dEcouDH[resid(i, resLig, residues)]))
            write(resMM_VDW_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEvdw[resid(i, resLig, residues)], dEvdw[resid(i, resLig, residues)]))
            write(resMM_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dEcou[resid(i, resLig, residues)] + dEvdw[resid(i, resLig, residues)], dEcouDH[resid(i, resLig, residues)] + dEvdw[resid(i, resLig, residues)]))
            write(resPBSA_pdb, @sprintf("%s%6.2f%6.2f\n", atomline, dPBres[resid(i, resLig, residues)] + dSAres[resid(i, resLig, residues)], dPBres[resid(i, resLig, residues)] + dSAres[resid(i, resLig, residues)]))
        end

        close(resMM_pdb)
        close(resMM_COU_pdb)
        close(resMM_VDW_pdb)
        close(resPBSA_pdb)
        close(resPBSA_PB_pdb)
        close(resPBSA_SA_pdb)
        close(res_MMPBSA_pdb)

        t_end = now()
        remain_time = (Nfrm - fr) * (t_end - t_start).value / 1000
        @printf("  MM-PBSA(with DH) = %9.3f(%9.3f) kJ/mol    Estimated Time Remaining: ~ %0dh:%02dm:%02ds\n",
                    dh[fr, 1], dh[fr, 2], Int32(floor(remain_time / 3600)), Int32(floor(remain_time % 3600 / 60)), remain_time % 60)
    end

    res_dat = open(joinpath(tprpath, "$pid~res.dat"), "w")
    write(res_dat, "#mol res       MM-PBSA(with DH  )        MM(with DH  )      PBSA       COU(with DH  )       VDW        PB        SA\n")
    for i = 1:Nres
        VDWres[i]    /= Nfrm
        PBres[i]     /= Nfrm
        SAres[i]     /= Nfrm
        COUres[i, 1] /= Nfrm
        COUres[i, 2] /= Nfrm

        MMres[i, 1] = COUres[i, 1] + VDWres[i]
        MMres[i, 2] = COUres[i, 2] + VDWres[i]
        dHres[i, 1] = MMres[i, 1] + PBres[i] + SAres[i]
        dHres[i, 2] = MMres[i, 2] + PBres[i] + SAres[i]

        resname = residues[i]
        resname = replace(resname, r"L~0+" => "Lig ")
        resname = replace(resname, r"P~0+" => "Pro ")
        write(res_dat, @sprintf("%-12s %9.3f(%9.3f) %9.3f(%9.3f) %9.3f %9.3f(%9.3f) %9.3f %9.3f %9.3f\n", 
        resname, dHres[i, 1], dHres[i, 2], MMres[i, 1], MMres[i, 2], PBres[i] + SAres[i], COUres[i, 1], COUres[i, 2], VDWres[i], PBres[i], SAres[i]))

        ii = ""
        if i == 1
            ii = @sprintf("---------------------------------------\n%-9s", "mean")
        end
        endres = ""
        if i == Nres
            endres="\n"
        end
        if withLig == 1
            write(res_MMPBSA, @sprintf("%s%9.3f(%9.3f) %s", ii, dHres[i, 1], dHres[i, 2], endres))
            write(resMM, @sprintf("%s%9.3f(%9.3f) %s", ii, MMres[i, 1], MMres[i, 2], endres))
            write(resMM_VDW, @sprintf("%s%9.3f %s", ii, VDWres[i], endres))
            write(resMM_COU, @sprintf("%s%9.3f(%9.3f) %s", ii, COUres[i, 1], COUres[i, 2], endres))
        end
        write(resPBSA, @sprintf("%s%9.3f %s", ii, preK*(PBres[i]+SAres[i]), endres))
        write(resPBSA_PB, @sprintf("%s%9.3f %s", ii, preK*(PBres[i]), endres))
        write(resPBSA_SA, @sprintf("%s%9.3f %s", ii, preK*(SAres[i]), endres))
    end

    close(res_dat)
        
    for i = 1:Npro
        ii = Ipro + i - 1
        atomline = @sprintf("ATOM  %5d %-4s %3s A%4d    %8.3f%8.3f%8.3f", ii, type[ii], res[ii], resID[ii], x[ii], y[ii], z[ii])
        if withLig == 1
            write(res_MMPBSA, @sprintf("%s%6.2f%6.2f\n", atomline, dHres[resid(i, resPro, residues), 1], dHres[resid(i, resPro, residues), 2]))
            write(resMM, @sprintf("%s%6.2f%6.2f\n", atomline, MMres[resid(i, resPro, residues), 1] + VDWres[resid(i, resPro, residues)], MMres[resid(i, resPro, residues), 2] + VDWres[resid(i, resPro, residues)]))
            write(resMM_COU, @sprintf("%s%6.2f%6.2f\n", atomline, COUres[resid(i, resPro, residues), 1], COUres[resid(i, resPro, residues), 2]))
            write(resMM_VDW, @sprintf("%s%6.2f%6.2f\n", atomline, VDWres[resid(i, resPro, residues)], VDWres[resid(i, resPro, residues)]))
        end
        write(resPBSA, @sprintf("%s%6.2f%6.2f\n", atomline, preK * (PBres[resid(i, resPro, residues)] + SAres[resid(i, resPro, residues)]), preK * (PBres[resid(i, resPro, residues)] + SAres[resid(i, resPro, residues)])))
        write(resPBSA_PB, @sprintf("%s%6.2f%6.2f\n", atomline, preK * PBres[resid(i, resPro, residues)] + SAres[resid(i, resPro, residues)], preK * PBres[resid(i, resPro, residues)] + SAres[resid(i, resPro, residues)]))
        write(resPBSA_SA, @sprintf("%s%6.2f%6.2f\n", atomline, preK * SAres[resid(i, resPro, residues)], preK * SAres[resid(i, resPro, residues)]))
    end
    for i = 1:Nlig
        ii = Ilig + i - 1
        atomline = @sprintf("ATOM  %5d %-4s %3s A%4d    %8.3f%8.3f%8.3f", ii, type[ii], res[ii], resID[ii], x[ii], y[ii], z[ii])
        write(res_MMPBSA, @sprintf("%s%6.2f%6.2f\n", atomline, dHres[resid(i, resLig, residues), 1], dHres[resid(i, resLig, residues), 2]))
        write(resMM, @sprintf("%s%6.2f%6.2f\n", atomline, MMres[resid(i, resLig, residues), 1] + VDWres[resid(i, resLig, residues)], MMres[resid(i, resLig, residues), 2] + VDWres[resid(i, resLig, residues)]))
        write(resMM_COU, @sprintf("%s%6.2f%6.2f\n", atomline, COUres[resid(i, resLig, residues), 1], COUres[resid(i, resLig, residues), 2]))
        write(resMM_VDW, @sprintf("%s%6.2f%6.2f\n", atomline, VDWres[resid(i, resLig, residues)], VDWres[resid(i, resLig, residues)]))
        write(resPBSA, @sprintf("%s%6.2f%6.2f\n", atomline, preK * (PBres[resid(i, resLig, residues)] + SAres[resid(i, resLig, residues)]), preK * (PBres[resid(i, resLig, residues)] + SAres[resid(i, resLig, residues)])))
        write(resPBSA_PB, @sprintf("%s%6.2f%6.2f\n", atomline, preK * PBres[resid(i, resLig, residues)] + SAres[resid(i, resLig, residues)], preK * PBres[resid(i, resLig, residues)] + SAres[resid(i, resLig, residues)]))
        write(resPBSA_SA, @sprintf("%s%6.2f%6.2f\n", atomline, preK * SAres[resid(i, resLig, residues)], preK*SAres[resid(i, resLig, residues)]))
    end

    close(resPBSA)
    close(resPBSA_PB)
    close(resPBSA_SA)
    close(resMM_VDW)
    close(resMM)
    close(resMM_COU)
    close(res_MMPBSA)

    dH = [0.0, 0.0]
    MM = [0.0, 0.0]
    COU = [0.0, 0.0]
    VDW = 0
    PB = 0
    SA = 0
    for i = 1:Nfrm
        dH[1] += dh[i, 1] / Nfrm
        dH[2] += dh[i, 2] / Nfrm
        MM[1] += mm[i, 1] / Nfrm
        MM[2] += mm[i, 2] / Nfrm
        COU[1] += cou[i, 1] / Nfrm
        COU[2] += cou[i, 2] / Nfrm
        VDW += vdw[i] / Nfrm
        PB += pb[i] / Nfrm
        SA += sa[i] / Nfrm
    end

    TdS = [0.0, 0.0]
    dG = [0.0, 0.0]
    Ki = [0.0, 0.0]
    if withLig == 1
        for i = 1:Nfrm
            TdS[1] += exp((mm[i, 1] - MM[1]) / RT2kJ) / Nfrm
            TdS[2] += exp((mm[i, 2] - MM[2]) / RT2kJ) / Nfrm
        end
    end
    TdS[1] = -RT2kJ * log(TdS[1])
    TdS[2] = -RT2kJ * log(TdS[2])
    dG[1] = dH[1] - TdS[1]
    dG[2] = dH[2] - TdS[2]
    Ki[1] = exp(dG[1] / RT2kJ)
    Ki[2] = exp(dG[2] / RT2kJ)

    write(MMPBSA, @sprintf("----------------------------------|------------------------------------------|--------------------------------|
%-12s %9.3f(%9.3f) | %9.3f(%9.3f) %9.3f %9.3f | %9.3f(%9.3f) %9.3f |
%-12s %9.3f(%9.3f) |
----------------------------------|
%-12s %9.3f(%9.3f) kJ/mol = %9.3f(%9.3f) kcal/mol
%-12s %9.3E(%9.3E)  uM    = %9.3E(%9.3E)    nM\n", 
"  dH=",   dH[1],     dH[2],     MM[1], MM[2], PB, SA, COU[1], COU[2], VDW,
"-TdS=", -TdS[1],   -TdS[2],
"  dG=",   dG[1],     dG[2],     dG[1]/4.184, dG[2]/4.184,
"  Ki=",   Ki[1]*1E6, Ki[2]*1E6, Ki[1]*1E9,   Ki[2]*1E9))

    write(MMPBSA, "+==============================================================+
   dG = RTln(Ki) = -RTln(KA)                          
   Ki = 1/KA = exp(dG/RT) = IC50 = EC50               
+==============================================================+
| M(mol/L) |  m/u/pM | dG(kcal/mol) | dG(kJ/mol) | Affinity    |
|==============================================================|
| 0.1      | 100  mM |    -1.364    |   -5.708   |             |
| 0.01     |  10  mM |    -2.728    |  -11.416   |             |
| 0.001    |   1  mM |    -4.093    |  -17.124   |             |
| 0.0001   | 100  uM |    -5.457    |  -22.832   | Weak        |
|--------------------------------------------------------------|
| 0.00001  |  10  uM |    -6.821    |  -28.540   |             |
| 1.0E-06  |   1  uM |    -8.185    |  -34.248   | Medium      |
|--------------------------------------------------------------|
| 1.0E-07  | 100  nM |    -9.550    |  -39.956   |             |
| 1.0E-08  |  10  nM |   -10.914    |  -45.664   |             |
| 1.0E-09  |   1  nM |   -12.278    |  -51.372   | Strong      |
|--------------------------------------------------------------|
| 1.0E-10  | 100  pM |   -13.642    |  -57.080   |             |
| 1.0E-11  |  10  pM |   -15.007    |  -62.788   |             |
| 1.0E-12  |   1  pM |   -16.371    |  -68.496   | Very strong |
+==============================================================+\n")
    close(MMPBSA)

    finish_hint("Finished MMPBSA calculation.\n")
    rm(joinpath(wd, "io.mc"))
    rm(joinpath(tprpath, "mdout.mdp"))
    rm(joinpath(tprpath, "mdout.err"))
    # request_hint("Remove intermediate files? [Y/n]")
    # res = get_input("y")
    # if res == "Y" || res == "y"
    #     [rm(joinpath(wd, f)) for f in readdir(wd) if startswith(f, "$pid") && (endswith(f, "pqr") && endswith(f, "err"))]
    # end
end

end # module
